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Abstract 

These lectures given to graduate students in theoretical particle physics, provide an 
introduction to the "inner workings" of computer algebra systems. Computer algebra 
has become an indispensable tool for precision calculations in particle physics. A 
good knowledge of the basics of computer algebra systems allows one to exploit these 
systems more efficiently. 
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1 Introduction 

Many of the outstanding calculations leading to precise predictions in elementary particle 
physics would not have been possible without the help of computer algebra systems (CAS). 
Examples of such calculations are the calculation of the anomalous magnetic moment of 
the electron up to three loops [1], the calculation of the ratio 

„ e+e~ — > hadrons , , 



of the total cross section for hadron production to the total cross section for the production 
of a iJ.~^n~ pair in electron-positron annihilation to order O (a^) [2] (also involving a three 
loop calculation) or the calculation of the 4-loop contribution to the QCD /9-function [3] . 
While the examples cited above have been carried out to an impressive third or fourth or- 
der, they only involve a single scale. Additional complications occur if multiple scales are 
involved. This is a challenging problem and computational quantum field theory witnessed 
a tremendous boost in the past two years in this area. The breakthrough was the devel- 
opment of systematic algorithms, which allowed the calculation of higher loop amplitudes 
with more than one scale. Examples are here the calculation of the most interesting 2^2 
processes [4, 5, 6, 7], which depend on two scales and the calculation for e+e" —>■ 3 jets 
[8, 9], which depends on three scales. None of the examples cited above could have been 
done without the help of computer algebra. Turning the argument around, these calcu- 
lations are just the ones which are at the edge of what is feasible with todays methods 
and technology. Limiting factors are not only missing knowledge how to calculate certain 
processes, but also the speed and memory of computers, missing knowledge on efficient 
algorithms or improper representation of the data inside the computer. Understanding the 
basics of computer algebra systems allows one to exploit these systems more efficiently. 
This is the main goal of these lectures. 



Particle physics is one important field of application for computer algebra and exploits 
the capabilities of computer algebra systems. This leads to valuable feed-back for the 
development of computer algebra systems. Quite a few computer algebra systems have 
their roots within the high energy physics community or strong links with them: RE- 
DUCE [10], SCHOONSHIP [11], MATHEMATICA [12], FORM [13] or GiNaC [14], to 
name only a few. Looking at the history of computer algebra systems, the first programs 
date back to the 1960's. Fig. (1) gives a historical overview together with the first ap- 
pearance of some programming languages. The first systems were almost entirely based 
on LISP ("List Programming language"). LISP is an interpreted language and, as the 
name already indicates, designed for the manipulation of lists. Its importance for symbolic 
computer programs in the early days can be compared to the importance of FORTRAN 
for numerical programs in the same period. Already in this first period, the program RE- 
DUCE had some special features for the application to high energy physics. An exception 
to the LISP-based programs was SCHOONSHIP, written in assembler language by Velt- 
man and specially designed for applications in particle physics. The use of assembler code 
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lead to an incredible fast program (compared to the interpreted programs at that time) 
and allowed the calculation of more complex scattering processes in high energy physics. 
The importance of this program was recognized in 1998 by awarding the Nobel prize to M. 
Veltman. Also the program MACSYMA [15] deserves to be mentioned explicitly, since it 
triggered important development with regard to algorithms. In the 1980's new computer 
algebra systems started to be written in C. This allowed to exploit better the resources 
of the computer (compared to the interpreted language LISP) and at the same time al- 
lowed to maintain portability (which would not have been possible in assembler language). 
This period marked also the appearence of the first commercial computer algebra system, 
among which Mathematica and Maple [16] are the best known examples. In addition, also 
a few dedicated programs appeared, an example relevant to particle physics is the program 
FORM by J. Vermaseren as a (portable) successor to SCHOONSHIP. In the last few years 
issues of the maintainability of large projects became more and more important and the 
overall programming paradigma changed from procedural programming to object-oriented 
design. In terms of programming languages this was reflected by a move from C to C-I--I-. 
Following this change of paradigma, the library GiNaC was developed. The GiNac library 
allows symbolic calculations in C++. At the same time the last few years marked also a 
transition away from commercial systems (back) to open-source projects. 

The choice of the approriate computer algebra system depends on the personal needs. 
It is worth analyzing the problem to be solved first. One class of problems requires only 
basic support from the computer algebra systems. These problems are often "local" prob- 
lems, where the problem expands into a sum of different terms and each term can be solved 
independently of the others. The complication which occurs is given by the fact that the 
number of terms can become quite large. Here the requirements on a computer algebra 
system are bookkeeping and the ability to handle large amounts of data. Quite a few 
problems in high energy physics fall into this class. A notorious example are tree-level 
processes involving several triple gauge-boson vertices. The Feynman rules, like the one 
for the three-gluon vertex, 

= gr'' [g"' {k's - k^2) + g'' iK - k',) + g'" {k^ - k',)] , (2) 

expand this vertex into up to six terms, which easily blow up the size of intermediate 
expressions. 

The second class of problems requires more sophisticated methods. These can be either 
standarized non-local operations (factorization is an example), which to some extend are 
already implemented in some computer algebra systems, or dedicated algorithms, which are 
developed by the user to solve this particular problem. Here the ability to model abstract 
mathematical concepts in the programming language of the computer algebra system is 
essential. 

Certainly, the two complications can also occur at the same time, e.g. the need to im- 
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Figure 1: A historical survey on computer algebra system together with the first appearance 
of some programming languages. 
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plement dedicated algorithms and the need to process large amounts of data. Higher 
loop calculations in quantum field theory tend to fall into this category. To summarize, 
requirements on a computer algebra system are therefore: 

• Efficiency in performance: If the system needs to process large amounts of data, 
performance is a priority. Usually this implies that the user's program is compiled 
and not interpreted. Systems optimized for performance may also contain low-level 
routines, which exploit efficiently the resources of the computer (memory, hard disc). 

• Support of object-oriented programming: This is of importance for users, who would 
like to implement complex algorithms for abstract mathematical entities. Being able 
to program at an abstract level (e.g. not in a low-level computer language) can 
reduce significantly the number of bugs. In a similar direction goes the requirement 
for strong type checking. Strong type checking can catch a number of bugs already 
at compilation time. 

• Short development cycles: It is usually the case that most time is spent for the 
development and the implementation of algorithms. The actual running time of 
the completed program is usually negligible against the development time. Since 

computer programs are developed by humans, the development time can be reduced 
if the computer algebra system allows for interactive use and provides high-quality 
output on the screen. 

• Source code freely available: No computer program is free of bugs, neither commer- 
cial products nor non-commercial programs. The probability of finding a bug in a 
particular computer algebra system while working on a large project with this com- 
puter algebra system is not negligible. If the source code is available it is simpler to 
trace down the bug, understand its implications and to fix it. 

Unfortunately, there is no system which would fuUfiU all requirements and the choice of the 
appropriate system involves therefore some inevitable trade-off. As examples for computer 
algebra systems I summarize the main features of Maple, FORM and GiNaC: 

• Maple: This is a commercial product. The advantages are a graphical user interface 
and the support for some sophisticated operations like factorization, integration or 
summation. The disadvantages are that Maple is quite inappropriate when it comes 
to processing large amounts of data. 

• FORM: The advantage of FORM is its speed and its capability to handle large 
amounts of data. It is widely used within the high energy community and the program 
of choice for calculations involving large intermediate expressions. 

• GiNaC: This is a library for symbolic computations in C++. The main feature is 
the possibility to implement one's own algorithms in an object-oriented way. It can 
handle large amounts of data and for some benchmark tests the speed is comparable 
to FORM. 
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These lecture notes follow a bottom-up approach from basic data structures to complex 
algorithms for loop calculations in high energy physics. In the next section I discuss data 
structures, object-oriented programming and the design of a simple toy computer algebra 
system. The discussion of the design will follow closely the design of the GiNaC-library. 
This is partly motivated by the fact that here the source code is freely available, as well as 
partly by the fact that I am involved in the development of a program [17], which combines 
the GiNaC-library with a C/C++ interpreter and a what-you-see-is-what-you-get editor. 
The resulting system allows an interactive use and displays results in met afont- type qual- 
ity. 

Section 3 deals entirely with efficiency. The efficient implementation of algorithms is dis- 
cussed in detail for three examples: Shuffling relations, the multiplication of large numbers 
and the calculation of the greatest common denominator. 

The development of computer algebra systems initiated also research on systematic algo- 
rithms for the solution of certain mathematical problems. In section 4 I discuss some of 
the most prominent ones: Factorization, symbolic integration, symbolic summation and 
simplifications with the help of Grobner bases. 

Section 5 is devoted to computational perturbative quantum field theory. I review algo- 
rithms, which have been developed recently, for the calculation of multi-loop integrals. 
General textbooks, on which parts of these lectures are based, are: D. Knuth, "The Art of 
Computer Programming" [18], K. Ceddcs ct al., "Algorithms for Computer Algebra" [19] 
and J. von zur Gathen and J. Gerhard, "Modern Computer Algebra" [20]. 

2 Data structures 

In this section I start from the basics: How data structures representing mathematical 
expressions, can be stored in the physical memory of the computer. Additional information 
on this subject can be found in the excellent book by Knuth [18]. The first programs for 
symbolic computations implemented mathematical expressions by nested lists and I discuss 
as an example symbolic differentiation in LISP first. Lists are special cases of "containers", 
and alternative types of containers are discussed in the following. In particular I examine 
the CPU time necessary to access or operate on the elements in various containers. The size 
of programs for the solution of certain problems depends on the complexity of the problem 
and can become large. In order to avoid that these programs become unmanageable, 
object-oriented techniques have been invented. I review these techniques briefiy. At the 
end of this section I discuss the design of a simple toy computer algebra system. 

2.1 Lists 

I start the section on data structures with a concrete example: Symbolic differentiation 
programmed in LISP. This is THE classic example for symbohc calculations. Symbolic 
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differentiation can be specified by a few rules: 

d 

—c = 0, 
ax 

—X = 1, 
ax 

^(/WsW) = + (3) 

These rules are sufficient to differentiate polynomials in x. An implementation in LISP^ 
looks like this: 

(DEFUN OPERATOR (LIST) (CAR LIST)) 



(DEFUN ARGl (LIST) (CADR LIST)) 
(DEFUN ARG2 (LIST) (CADDR LIST)) 



(DEFUN DIFF (E X) 

(COND ((ATOM E) (COND ((EQUAL E X) 1) 

(T 0))) 
((EQUAL (OPERATOR E) '+) 

'(+ ,(DIFF (ARGl E) X) , (DIFF (ARG2 E) X))) 
((EQUAL (OPERATOR E) '*) 
'(+ (* ,(DIFF (ARGl E) X) , (ARG2 E)) 

(* ,(ARG1 E) ,(DIFF (ARG2 E) X)))))) 

A few comments on this short program are in order. Within LISP the prefix notation is 
usually used, e.g. A + B is represented by (+ A B) . The first three lines define aliases to the 
build-in LISP functions CAR, CADR and CADDR to extract the first, second or third element 
of a list, respectively. They are only introduced to make the program more readable. 
The program itself shoud be readable even without the knowledge of LISP, except for the 
appearance of the single quote " ' the backquote " ' " and the comma "," character. 
LISP generally interprets the first element in a list as the name of an operation and tries 
to evaluate this operation with the remaining elements of the list as arguments. This 
behaviour can be prohibited by putting a single quote " ' " in front of the list and the 
list remains unevaluated. The backquote " ' " acts like the single quote, except that any 
commas that appear within the scope of the backquote have the effect of unquoting the 
following expression. 

In this simple example it is further assumed that addition and multiplication take only two 
arguments. LISP is an interactive language and entering 

^For the nostalgic: A free LISP interpreter is available from http://www.gnu.org/software/gcl/gcl.htnil. 
A good introduction to LISP and programming techniques in LISP can be found in [21]. 
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(DIFF A X) ^X) 



at the prompt for 




(4) 



yields 



(+ (* OX) (* A D), 



which stands for 



(O-x) + (a - 1). 



(5) 



Simphfications hke 0-a; = 0, a-1 = a or 0+a = a are out of the scope of this simple example. 
This example already shows a few important features of computer algebra programs: 

• The distinction between objects, which contain sub-objects and objects without fur- 
ther substructures. The former are generally referred to as "containers", the later 
are called "atoms". Examples for atoms are symbols like a or x, while examples 
for containers are data structures which represent multiplication or addition of some 
arguments. 

• The use of recursive techniques. In the example above the function DIFF calls itself 
whenever it encounters a multiplication or addition. Note that the recursive function 
call has simpler arguments, so that the recursion will terminate. 

• Lists are used to represent data structures hke addition and multiplication. The first 
element in the list specifies what the list represents. 

• Lists can be nested. The output of the example above, (+ (* X) (* A 1)), con- 
sists of two hsts nested inside another hst. 

• The output of a symbolic procedure is not necessarily in the most compact form. 
2.2 Containers 

A container is an object that holds other objects. Lists and arrays are examples of contain- 
ers. There are several ways how the information of a container can be stored in physical 
memory. In particular, the time needed to access one specific element will depend on the 
lay-out of the data in the memory. The appropriate choice depends on the specific problem 
under consideration. In this context, the "big-0" notation is useful: An indication 0{n) 
for an operation on n elements means that the operation takes time proportional to the 
number of elements involved. Fig. (2) gives a rule-of-thumb for the significance of the cost 
of certain operations. 

An array (also called "vector" within the C-|— I- terminology) stores the information hnearly 
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0(1) 

o(iogH) 

0{n) 

0{n log(n)) 



cheap 

fairly cheap 
expensive 

expensive 
very expensive 



Figure 2: A summary on the costs of certain operations. Operations of order 0(1) or 
0(log(n)) are considered "cheap" operations. Further, 0{n\og{n)) is considered closer to 
0{n) than to 0{n^). It is usually a considerable speep-up, if an operation which naively 
takes O(n^) time, can be improved to 0(nlog(n)). In generally one tries to avoid opera- 
tions, which take O(n^) time. 

in memory. Given the fixed size / of one entry and the address Oi of the first entry, the 
address aj of the i-th entry is easily obtained as 



This involves one multiphcation and one addition. The CPU time to access one specific 
element is therefore independent of the number of elements in the array. However, suppose 
that we have an array of n elements and that we would like to insert a new element between 
the ith and {i + l)th entry. This involves to shift the entries with index n, n — 1, n — 2, 
i + 1 by one position and is an 0(n) operation. If these operations occur frequently, an 
array is not the best suited data structure. 

In that case a list structure is more appropriate. A list is often implemented as a double- 
linked list, where each node of the double-linked list contains one field of information, one 
pointer to the next element and one pointer to the previous element. Inserting a new 
element somewhere in the middle of the list involves only updating the pointers and is an 
0(1) operation. Typical fist operations are to add or to remove elements in the middle of 
the list. However there is also a drawback for lists. To access the ith element of a list, 
there is no other way than to start at the first element and to follow the pointers to the 
next elements sequentially, until one arrives at the ith entry. This is an 0(n) operation. 
A generalization of a double-linked list is a rooted tree, where each node may have several 
sub-nodes. A list can be viewed as a tree, where each node has exactly one sub-node. 
An important special case is a binary tree, where each node has up to two sub-nodes. A 
binary tree can be used to encode an order relation: Elements of the left sub-tree of a 
particular node are "less-than" the current node, whereas elements of the right sub-tree 
are "greater-than" the current node. 

An associative array (sometimes also called a map) keeps pairs of values. Given one value, 
called the key, one can access the other, called the mapped value. One can think of an 
associative array as an array for which the index need not be an integer. An associative 
array can be implemented by a binary tree. In this case it is further assumed that there 
is a less-than operation for the keys and the associative array keeps it's elements ordered 
with respect to this relation. To find the mapped value corresponding to a specific key k, 



ai + (i - 1)1. 



(6) 
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Figure 3: The cost for specific operations with different types of containers. A dash indicates 
that the corresponding operation is usually not provided for this container. For the hash 
map it is assumed that the map is not tightly filled, such that collisions occur rarely. 



one first compares the key k with the key kn of the entry at the root of the tree. If A; < kn, 
the mapped value is in the left sub-tree and the procedure is repeated with the top-node of 
the left sub-tree as root, if k > kr the mapped value is in the right sub-tree and \i k = k^ 
we have already found the corresponding (key, value) pair. To find the mapped value takes 
on the average 0(log(n)) operations. It is important to note that insertion of new elements 
is also an 0(log(n)) operation and does not require to sort the complete tree. Insertion is 
done as follows: One starts at the root and one compares the key k of the new clement 
with the key k^ of the root. Let us assume that k > k^. If the root has a right sub-tree, 
the procedure is repeated with the top-node of the right sub-tree as root. If the root does 
not have a right sub-tree, we attach a new right sub-tree to the root, consisting of the new 
element. 

An alternative implementation for an associative array uses a hash map. Suppose that 
there is an easy to calculate function, called hash function, which maps each key to a spe- 
cific address in a reserved memory area. In general, there will be no hash function, which 
can guarantee that two different keys are not mapped to the same address, but a good 
hash function will avoid these coUisions as much as possible. If a collision occurs, there 
are several strategies to deal with this case: The simplest implementation for insertions 
would just use the next free entry, whereas to access elements one would do a linear search. 
An associative array based on a hash map requires therefore a hash function and an "is 
equal" -operation for its keys. Since the calculation of the hash value should be a fast oper- 
ation, many hash functions are based on exclusive-or operations and bit-wise rotations. If 
collisions are not frequent, insertions and access can be done in constant time. Collisions 
arc less frequent, if the memory area for the hash map is not tightly filled. Therefore a hash 
map is appropriate if speed for access and insertions is important and sufficient memory is 
available. 

Fig. (3) summarizes the cost for specific operations involving the various types of contain- 
ers. Back operations are operations, where access, insertion or deletion occur at the end 
of the container. A typical example are stack operations. 
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2.3 Object-oriented design: A little bit C-\ — |- 

Objcct-oricntcd programming techniques were invented to allow the development and main- 
tenance of large, complex software applications. Experience showed that techniques which 
worked for smaller projects do not necessarily scale to larger projects. In this section I dis- 
cuss some features of C-|— I- as an example for an object-oriented programming language. A 
detailed introduction to C++ and object-oriented techniques can be found in the book by 
Stroustrup [22]. Most of the topics discussed here are also present in other object-oriented 
languages or can be modelled by the user in languages which do not have native support 
for object-oriented techniques. 

There are several reasons, which motivate the particular choice of C++ from the set of 
object-oriented programming languages: Since 1998 C++ is standardized, which enhances 
the portability of programs to different platforms. This is of particular importance in 
academia, where one moves in the early stages of a career from one post-doc position to 
another. Furthermore, C++ is widely used. As a consequence, additional development 
tools are available. The hst includes debuggers, editors with automatic high-lighting facih- 
ties for C++ code and tools for the automatic generation of documentation. Finally, C++ 
allows operator overloading. In the scientific domain, operations like addition and multi- 
plication are frequently used. Operator overloading allows one to use the same notation 
for different data types, e.g. one can write for example 

2*a + b 

independently if a and b are numbers, vectors or matrices. This makes programs more 
readable. It is also a major argument against Java in scientific computing. Java does not 
allow operator overloading and writing 

V . add_vec (w) ; 

for the addition of two vectors makes programs less readable. 

A first tool for object-oriented programming is a modular approach corresponding to a "di- 
vide and conquer" strategy. If a problem can be clcary separated into two sub-problems, 
the complexity is already greatly reduced. A basic module or the smallest independent 
entity in C++ is a class. A class consists of data members and methods operating on 
the data (methods are also called member functions). An example for a class would be 
the complex numbers, whose data members are two double precision variables x and y, 
representing the real and imaginary part, respectively. Member functions could be a print 
routine, addition and subtraction, etc.. Within the "divide and conquer" strategy falls 
also the strict separation of the implementation from the interface. This is done by private 
and public members. Data members are usually taken to be private members, to protect 
them to be changed accidently from outside the class. A well designed separation between 
an interface and an implementation allows one to replace the implementation of a specific 
method by a more efficient algorithm without changing the interface. Therefore no changes 
have to be done to the rest of the program. 
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A second tool for object-oriented programming is data abstraction. General purpose pro- 
gramming languages come usually with a limited number of build-in data types (integers, 
real double precision, etc.). In many cases the user would like to have additional data 
types, in high energy physics one might like to have a data type "complex number" or 
even a data type "Feynman diagram" . Now, there might be programming languages which 
have a build-in data type "complex number", but it is quite unlikely to find a language 
which has a data type "Feynman diagram". The ability to extend the build-in data types 
by user-defined data types is therefore essential. Furthermore it occurs quite often that 
several data types arc related to each other. In C-I--I- similar concepts are modeled through 
inheritance and derived classes. For example, if there are the data types "strongly inter- 
acting particle", "gluon" and "quark", the type "strongly interacting particle" would be 
a base class from which the classes "gluon" and "quark" are derived. Derived classes are 
specializations of base classes and may have additional properties. Therefore, a derived 
class may be used where-ever only the base class is required. Inheritance and derived 
classes can be used to avoid that similar code is implemented more than once. Within the 
C-|--|- community one draws in inheritance diagrams arrows from the derived class pointing 
to the base class. 

A third tool of object-oriented programming is the ability to re-use and extend a given 
program. The traditional re-use of an existing program consists in an application, where 
the new code calls the old code. Object-oriented programming techniques allow also the 
reverse situation, e.g. that the old code calls the new code. In C-I--I- this can be done 
through virtual functions which are resolved at run-time. A typical example is the situ- 
ation where the old program is a larger framework, which one would like to extend with 
some new functionality. If implemented properly, this will work even without recompiling 
the old code. 

The fourth tool of object-oriented programming techniques are generic algorithms. For ex- 
ample, sorting a list of integers or a list of real double precision variables can be done with 
the same algorithm. In C-|— I- the use of templates allows the implementation of generic 
algorithms. The Standard Template Library (STL) already contains implementations for 
the most important algorithms, among other things it provides vectors, lists, maps and 
hash maps, discussed in the previous section. 

C-I--I- was developed from the C programming language and retains C as a subset. As in 
C, data can be stored in memory either automatically on the stack, dynamically allocated 
on the heap or statically at a fixed address. For subroutine calls it is usually not efficient 
to copy a large data structure to local variables upon entering the subroutine. Instead it 
is more appropriate to pass a reference or a pointer to the subroutine. 

2.4 A simple toy computer algebra system 

In this section I will discuss the design of a simple toy computer algebra system. This toy 
model should know about symbols ("a", "b", etc.), integers ("1", "2", etc.), addition and 
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multiplication. It should be possible to enter an expression like 

5a + 36 + 2a (7) 

which is automatically simplified to 7a + 36. The discussion of the implementation for this 
toy program follows closely the structure of the GiNaC-hbrary [14, 23]. For the GiNaC- 
library the source code is freely available and the interested reader is invited to study the 

techniques explained below in a real case example. 

From the description above we can identify five objects: symbols, integers, addition, mul- 
tiplication and expressions. Symbols and integers cannot have any subexpressions and are 
atoms, addition and multiplication are containers. An expression can be an atom or a 
container. The natural design is therefore to start from a base class "basic" , representing 
a basic expression, from which the classes "symbol", "numeric", "add" and "mul" are de- 
rived. The lay-out of the atomic classes is straightforward: The class "symbol" must store 
it's name (a string like "a" or "x"), whereas the class "numeric" must store it's numeric 
value. Since both are derived from the class "basic", the relevant part of code in C-|— |- 
could look as follows: 

class symbol : public basic 
{ 

protected: 

std:: string name; 

}; 

class numeric : public basic 
{ 

protected: 

int value; 

}; 

Before dealing with the lay-out of the remaining classes (in particular with the lay-out 
of the class "basic"), it is worth to consider memory management and efficiency first. 
Already in our simple toy example, expressions can be nested and can become large. It 
is inefficient to copy large expressions upon entering or leaving a subroutine. Passing 
a pointer is more efficient. However, programming with pointers at a high level is not 
elegant and error-prone. One therefore hides all pointer operations into an encapsulating 
class "ex": 

class ex 
{ 

public : 

basic *bp; 

}; 
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The user deals only with the class "ex" and never manipulates pointers directly. An 
instance of the class "ex" consists only of a pointer to a "basic" -object and is therefore 
extremely light-weight. When a new expression is created in a subroutine, this subroutine 
returns an "ex", e.g. a pointer to some "basic" -object. This requires that the "basic"- 
objcct is created dynamically on the heap. If it would be created on the stack, it will get 
automatically destroyed at the end of the subroutine, and the program will abort as soon 
as it tries to access this object outside the subroutine. However with dynamic memory 
allocation one has to ensure that objects, which are no longer needed, get deleted and that 
the occupied memory space is freed. Otherwise one would run out of memory soon. This 
can be done with a technique called "reference counting with copy-on- write semantics" [24]. 
The class "basic" has a counter, which keeps track of how many times this "basic" -object 
is pointed to: 

class basic 
{ 

private : 

unsigned ref count; 

}; 

It happens quite often that some expression is assigned to more than one variable. With 
reference counting, there is no need to store the same "basic" -object twice in memory. In 
the following lines of code 

ex el = 3*a + 2*b; 
ex e2 = el; 

both el and e2 point to the same "basic" -object and no copying takes place in the 
second hue of this code. After these lines have been executed, the ref count variable of 
the "basic" -object is equal to 2. However, copying is necesarry when one expression is 
changed: 

e2 = e2 + 4; 

This is called copy-on-write semantics. In this example e2 pointed first to the expression 
3a -|- 2b, and then to the new expression 3a -|- 26 -|- 4. The implementation ensures that the 
reference counter of 3a-|-26 is decreased by one, as soon as e2 points to the new expression 
3a + 26 + 4. Obviously, if the ref count variable of some "basic" -object equals zero, it is 
no longer referenced by any variable and therefore it can be safely deleted. This releases 
the memory occupied by the "basic" -object and avoids memory leaks. 

Let us now discuss the lay-out of the container classes "add" and "mul". As an exam- 
ple we consider again the expression 3a + 26 + 4. The naive representation would be a 
nested structure, where the top-level container is of type "add" with three summands. 
One summand is the numerical value 4, while the other two summands are of type "mul" 
and given by 3a and 26, respectively. Unfortunately this representation will result in 
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rather deep trees, which are slow to manipulate. Since products of the form "numerical 
coefficient" times "something else" occur quite often, it is advantageous to introduce an 
additional data type for this pair: 

class expair 
{ 

public : 

ex rest; // first member of pair, an arbitrary expression 
ex coeff ; // second member of pair, must be numeric 

}; 

These pairs are easy to manipulate, for example, if two pairs in a sum have the same value 
for the variable rest, one can simply add their numerical coefficients. The introduction 
of these pairs flattens expression trees signiflcantly. Our sum 3a + 26 + 4 can now be 
represented by a sequence of two pairs (a, 3), (6, 2) and the overall numerical constant 4. If 
one would introduce powers into the computer algebra system, a similar structure is seen 
for products, for example bx'^y'^z can be represented by the sequence of pairs {x,2), (y, 3) 
and (z, 1) and the overall numerical constant 5. Again, if two pairs in the product have the 
same value for the variable rest, one can simply add the numerical values of their variables 
coeff, which now represent the exponents. It is therefore appropriate to introduce for the 
containers "add" and "mul" a common base class "expairseq", which implements this 
data structure, e.g. a sequence of pairs together with an additional numerical coefficient: 

class expairseq : public basic 
{ 

protected: 

std: : vector<expair> seq; 
ex overall_coef f ; 

}; 

The sequence of pairs is implemented by using a vector from the Standard Template 
Library. The classes "add" and "mul" are then derived from "expairseq" : 

class add : public expairseq 

{ }; 



class mul : public expairseq 

{ }; 

Note that "add" and "mul" do not need any additional data members. Finally, we overload 
the operators + and *: 

const ex operator+ (const ex & Ih, const ex & rh) ; 
const ex operator* (const ex & Ih, const ex & rh) ; 
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Figure 4: The class structure of a simple toy computer algebra system. Long arrows point 
from derived classes to the base classes. A dashed arrow indicates that this class is a "smart 
pointer" to another class. 



These operators call the constructors of the classes "add" and "mul", respectively. Fig. 
(4) summarizes the class hierarchy of the simple toy computer algebra program. Note that 
the class "expair" is an internal helper class and not derived from "basic". 

We now discuss how to implement automatic simplifications. For example, we would like 
to have that 5a + 36 + 2a is automatically evaluated to 7a + 36. To this aim we introduce 
for the class "basic" a status flag "evaluated" and a method evalO: 

class basic 
{ 

public : 

virtual ex eval(void) const; 
protected: 

mutable unsigned flags; 



Every time a "basic" -object is assigned, this flag is checked. If the flag is not set, the 
method evalO is called. The base class "basic" has a trivial evaluation routine, which 
does nothing except setting the flag to evaluated. This implementation is also sufficient 
for the atomic classes "symbol" and "numeric". The method evalO is a virtual function 
and can be redefined in derived classes. Since all expressions are always accessed through 
pointers via the class ex the mechanism of virtual functions ensures that the appropriate 
method is called. In our example, the classes "add" and "mul" would redefine the virtual 
function eval: 

class add : public expairseq 



>; 
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{ 

public : 

ex eval(void) const; 

} 

An implementation of the method evalO for the class "add" would sort the elements of 
the sequence of "expair"s and combine element which have the same "rest". Note that 

automatic evaluation happens the first time an object is assigned, not the first time an 
object is constructed. This is due to the fact, that an object of type "add" like 2x — x 
evaluates to which is of type "symbol" and not of type "add" . 



3 Efficiency 

This section is devoted to issues related to efficiency. I discuss three examples: Shuffling 
relations, fast multiplication and the calculation of the greatest common divisor. Shuffling 
relations are discussed as an application of recursive techniques. For the multiplication 
of large numbers I outline three methods: the classical text-book method, Karatsuba's 
algorithm and a method based on fast Fourier transform. An algorithm for the calculation 
for the greatest common denominator is of central importance for many other algorithms. 
Therefore the Euclidean algorithm and improvements are discussed in detail. 



3.1 Shuffling 



Recursive techniques are often used in symbolic calculations. A classical example are the 
Fibonacci numbers defined by /(O) = 1, /(I) = 1 and 



f{n) 



/(n-l) + /(n-2) 



(8) 



for n > 2. Recursive procedures are easily implemented, but it should be noted that in 
most cases they provide not an efficient way to solve a problem. For example, the recursive 
definition of the Fibonacci numbers is certainly the most efficient way to create a table 
with the first n Fibonacci numbers, but a recursive approach is highly inefficient, if we just 
need one Fibonacci number /(n), if n is large. For the Fibonacci numbers a closed formula 
is known. 



n+1 



1 - 



n+1 



(9) 



and it is more efficient to use this closed formula in the latter case. 

With this warning, I will now discuss some issues related to the implementation of recursive 
procedures in computer algebra systems. As an example I will use shuffling relations. 
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Shuffling relations occur frequently in calculations. I discuss them for Euler-Zagier sums, 
which are defined by 



Zmu-,mk{n) — ^2 7^ ■ ■ ■ 7~^- (-'■^) 

n>ii>i2>--->ik>0 

Euler-Zagier sums form an algebra, that is to say that the product of two Euler-Zagier 
sums with the same upper summation limit n is again a finite sum of single Euler-Zagier 
sums. The product can be worked out recursively via the formula 

■^mi,...,mfc(^) ^ -^m'j^,...,mj ('^) 
n ^ 

^ X] 7^^"»2,...,mfc(«l - l)Zm[,...,m[{il " 1) 

' 1 

ii=l 1 

n ^ 

+ :;^"^'"i'-."»fc(^2 - l)^mi,...,m;(«2 - 1) 

22=1 *2 

n ^ 
1=1 

This multiplication can be implemented by a subroutine, which takes three lists res, argl 
and arg2 as arguments. When the routine is entered the first time, res is empty and 
argl and arg2 contain (mi, ...,mk) and (m'^, ...,m[), respectively. If argl (or arg2) is the 
empty list, we are basically done: We append the content of arg2 (or of argl) to res and 
return res. Otherwise we remove the first element from argl and append it to res and 
use recursion, this corresponds to the first line on the r.h.s of eq. (11). Similar, for the 
second (or third line) of eq. (11) one removes the first element from arg2, (or the first 
elements from argl and arg2) and appends it (or the sum of the two elements) to res. 
Most computer algebra systems provide a list data type within their framework. In GiNaC 
this type is called 1st and is a derived class from basic. It is tempting to implement the 
shuffle multiplication using this data type. In GiNaC this could look as follows: 

ex shuffle_mul (const ex & res, const ex & argl, const ex & arg2) ; 

where it is understood that res, argl and arg2 will always point to a 1st. However, such 
an approach can lead to a considerable performance loss, since the arguments res, argl 
and arg2 are now under the spell of the automatic evaluation procedure. The algorithm 
according to eq. (11) shuffles basically elements from argl and arg2 to res. At each step 
new lists are created with appended or removed elements, and the automatic evaluation 
procedure will check for each new list, if all of its elements are already evaluated. Usually all 
elements are already evaluated, therefore this procedure is a waste of computer resources. 
It is therefore better to use for the lists a temporary data structure which is not related 
to the automatic evaluation procedure. In our case, a vector containing the type ex would 
do the job: 
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typedef std: :vector<ex> exvector; 

ex shuffle_mul (const exvector & res, const exvector & argl, 

const exvector & arg2) ; 

Note that in the discussion of the GiNaC-hbrary we aheady crossed the data structure 
expair, which is not derived from the class basic and therefore not subject to the auto- 
matic evaluation procedure. 

Situations where this technique can be used occur quite frequently. Another example is 
the shuffle product of iterated integrals. Let us consider the integral 

J tl- Zi J t2- Z2 J tk- Zk 


A product of two such integrals can be reduced to a sum of single integrals according to 
G(zi,...,Zk;y) X G(zk+i,...,Zn;y) = ^ G{z^^i), z„(^n),y), (13) 

shuffles 

where the sum is over all permutations of ^i, 2;n, which keep the relative order of ^i, 
Zk and Zk+i, Zn fixed. Recursively, we can write for the product 

G{zi, Zk] y) X G{zk+u Zn] y) = (14) 

^ d ^ d 

/- G{Z2., Zk, t)G{Zk+l, Zn, t) + G{Zi, Zk, t)G{Zk+2, Zn, t) 
I — Zi J t — Zk+1 


and the algorithm can be implemented in complete analogy to the one for eq. (11). 
Similar considerations apply if we just want to transform one representation to another 
form. An examples is a change of basis from Euler-Zagier sums to harmonic sums, defined 

by 

The difference between Euler-Zagier sums and harmoinc sums is the upper summation limit 
for the subsums: i — 1 for Euler-Zagier sums and i for harmonic sums. The conversion 
between the two bases uses the formulae 

n ^ 

■^mi,...,mfe ('^) — ^ ^ • mi ■^rn2,...,irn,{}l) ~ ■^mi+m2,m3,...,mfe (''^) ) 
11=1 
n ^ 

Smi,...,mkiP') ~ ^y^j ■ mi ^m2,---,mui}l ~ 1) + Sm^j^m2,m3,...,mkiP') ■ (1^) 
11=1 

and are implemented by routines taking two arguments: 

ex Zsuin_to_Ssuin(const exvector & res, const exvector & arg) ; 
ex Ssuin_to_Zsuin(const exvector & res, const exvector & arg); 
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3.2 Multiplication of large numbers 



In this paragraph I discuss issues related to the efficiency for the multiplication of large 
numbers. Multiprecision arithmetic is needed in many applications. For the practitioner, 
there are several freely-available libraries which provide support for multiprecision num- 
bers. Examples are the libraries GMP [25], CLN [26] or NTL [27]. 

Let a be an integer with n digits in the base B, e.g. 

a = ao + aiB + ... + an-iB''''^. (17) 

For simplicity we assume that n is even (if n is odd we can just add a zero to the front). 
We can write a as 

a = UhB^'/'^ + ai, (18) 

where ah and ai have now n/2 digits in base B. Let b be another integer with n digits. 
The product ab can be calculated as 

ab = ahbhB'' + {ahk + aibh) B""/^ + aik (19) 

This is the classical method. To multiply two integers with n digits requires 4 multiplica- 
tions of integers with n/2 digits and four additions. Usually the computational cost of the 
additions can be neglected against the one for the multiplications. It is not too hard to 
see that the classical method is of order O(n^). 

The efficiency can be improved for large integers by rewritting eq. (19) as follows: 

ab = ahbhB''+[ahbh + aibi-iah-ai)ibh-bi)]B''/'' + aibi (20) 

This method requires only 3 multiplications of integers with n/2 digits and was invented 
by Karatsuba [28]. It can be shown that this algorithms grows hke 0{n^°^^^) ^ 0{n^-^^) 
with the number of digits n. 

There is even a faster algorithm developed by Schonhage and Strassen [29] and based 
on the fast Fourier transform. It is most easily explained for the multiplication of polyno- 
mials. A polynomial of degree n is usually represented by its coefficients af. 

a{x) = ao + aix + a2x'^ + ... + a-nx"^. (21) 

However, a polynomial of degree n is also uniquely defined by the values at n -|- 1 distinct 
points xq, xi, Xn- This representation is called the modular representation. Suppose 
that we would like to multiply a polynomial a{x) of degree n by a polynomial b{x) of degree 
m. Then the product is of degree n + m. Suppose that we know the polynomials a{x) and 
b{x) at n -|- m -|- 1 distinct points: 

a{x) = {a{xo),...,a{xn+m)) , 

b{x) = {b{xo),...,b{Xn+m))- (22) 
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Then 



( a{xo) ■ b{xo), a{xn+m) ■ b{^n+m) ) (23) 

defines uniquely the polynomial a{x)b{x) and requires only 0{n + m) operations. In the 
modular representation, multiplication of two polynomials of degree n is an 0{n) operation. 
However, this is only the cost if the objects are already in this particular representation. 
In general, we have to add the costs for converting the input polynomials to the modular 
representations and the cost for converting the result back to a standard representation. 
The method will be competitive only if these transformations can be done at a cost lower 
than, say, Karatsuba multiplication. Using a method based on the fast Fourier transform, 
it tiuns out that these transformations can be done efficiently and that the cost for these 
transformations is O(nlogn). Let us first consider the transformation from the standard 
representation 

a{x) = ao + oix + ... + a„_ia;""-^ (24) 

of a polynomial of degree less than n to the modular represenation. For simplicity we 
assume that n is even. We have the freedom to chose the n distinct evaluation points 
and a clever choice can reduce the required amount of calculation. A suitable choice is to 
evaluate the polynomial at the special points (1, a;, a;^, a;"~^}, where a; is a primitive 
n-th root of unity, e.g. 

u"^ = 1, but V 1 for < A; < n. (25) 

A few examples for primitive roots of unity: In C, i and —i are primitive 4-th roots of 
unity, —1 is not a primitive 4-th root of unity, since (—1)^ = 1. In Z17, 4 is a 4-th root of 
unity, since 4^ = 1 mod 17 and 4^ = 16 mod 17 as well as 4^ = 13 mod 17. 
The evaluation of the polynomial at this particular set of evaluation points is nothing than 
performing a discrete Fourier transform. Given a finite set of values (oq, ai, a„_i), the 
discrete Fourier transform maps those values to the set (oq, Oi, On-i)) where 

n-1 

= ^aju'^. (26) 
j=o 

An efficient way to calculate the values (ao, Oi, On-i) is given by the fast Fourier trans- 
form. Since we assumed that n is even, the n-th roots of unity satisfy 

= -u\ (27) 

We write the polynomial a{x) in the form 

a{x) = h{x^)+x-c{x'^), (28) 

where the polynomials h{y) and c(|/) have at most one-half the degree of a{x). Since we 
have (0;'+"/^)^ = {u^)^ it is sufficient to evaluate the polynomials h{y) and c{y) at n/2 
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distinct points instead of n. This is the basic building block for the fast Fourier transform. 
Using this technique recursively yields a method which performs the Fourier transform at 
O(nlogn). In complete analogy, the conversion from the modular representation to the 
standard representation is done by the inverse discrete Fourier transform. Given a finite 
set of values {aQ,ai, ...,a„_i), the inverse discrete Fourier transform maps those values to 
the set ( ), where 

^ n— 1 

aj = -J2^kOO~^''- (29) 

^ k=0 

Again, a similar decomposition like in eq. (28) is used to speed up the calculation. 



3.3 The greatest common divisor and the Euchdean algorithm 

It is often required to simplify rational functions by canceling common factors in the 
numerator and denominator. As an example let us consider 

ftt'rt = (--yf- (30) 

One factor of {x + y) is trivially removed, the remaining factors are cancelled once we 
noticed that {x^ — y^) = {x + y){x — y). For the implementation in a computer algebra 
system this is however not the way to proceed. The factorization of the numerator and 
the denominator into irreducible polynomials is a very expensive calculation and actually 
not required. To cancel the common factors in the numerator and in the denominator it 
is sufficient to calculate the greatest common divisor (gcd) of the two expressions. The 
efficient implementation of an algorithm for the calculation of the gcd is essential for many 
other algorithms. Like in the example above, most gcd calculations are done in polynomial 
rings. It is therefore useful to recall first some basic definitions from ring theory: 
A commutative ring {R, +, •) is a set R with two operations + and •, such that {R, +) is 
an abelian group and • is associative, distributive and commutative. In addition we always 
assume that there is a unit element for the multiplication. An example for a commutative 
ring would be Zg, e.g. the set of integers modulo 8. In this ring one has for example 
3 + 7 = 2 and 2-4 = 0. From the last equation one sees that it is possible to obtain zero 
by multiplying two non-zero elements. 

An integral domain is a commutative ring with the additional requirement 

a-b — =^ a = 0or6 = (no zero divisors). (31) 
Sometimes an integral domain D is defined by requiring 

a • b — a • c and a ^ =^ b — c (cancellation law). (32) 
It can be shown that these two requirements are equivalent. An example for an integral 
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Euclidean domain 
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Figure 5: Hierarchy of domains. Upward pointing arrows indicate that a former domain is 
automatically also a latter domain. 



domain would be the subset of the complex numbers defined by 

S = ^^a + biV5 a,6,Gz| (33) 

An element u E D is called unit or invertible if u has a multiplicative inverse in D. The 
only units in the example eq. (33) are 1 and —1. We further say that a divides b if there 
is an element x G -D such that b = ax. In that case one writes a\b. Two elements a,b E D 
are called associates if a divides b and b divides a. 

We can now define the greatest common divisor (gcd): An element c G -D is called the 
greatest common divisor of a and b if c\a and c\b and c is a multiple of every other element 
which divides both a and b. Closely related to the gcd is the least common multiple (1cm) 
of two elements a and b: d is called least common multiple of a and b if a\d and b\d and d 
is a divisor of every other element which is a multiple of both a and b. Since gcd and 1cm 
are related by 

1cm a, 6 = 34 
gcd(a, b) 

it is sufficient to focus on an algorithm for the calculation of the gcd. 



An element p E D — {0} is called prime if from p\ab it follows that p\a or p\b. An el- 
ement p E D — {0} is called irreducible if p is not a unit and whenever p = ab either a or 
6 is a unit. In an integral domain, any prime element is automatically also an irreducible 
element. However, the reverse is in general not true. This would require some additional 
properties in the ring. 

An integral domain D is called a unique factorization domain if for all a G -D — {0}, either 
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Figure 6: Structure of polynomial rings in one variable and several variables depending on 
the underlying coefficient ring R. 



a is a unit or else a can be expressed as a finite product of irreducible elements such that 
this factorization into irreducible elements is unique up to associates and reordering. It 
can be shown that in an unique factorization domain the notions of irreducible element 
and prime element are equivalent. In a unique factorization domain the gcd exists and 
is unique (up to associates and reordering). The integral domain S in eq. (33) is not a 
unique factorization domain, since for example 

21 = 3-7= (l- 2iV5^ (^1 + 2iVtj (35) 

are two factorizations into irreducible elements. An example for a unique factorization 
domains is the polynomial ring Z[x] in one variable with integer coefficients. 

An Euclidean domain is an integral domain D with a valuation map u:!) — {0}— >N into 
the nonnegative integer numbers, such that v{ab) > v{a) for all a, 6 e D — {0}, and for all 
a,b e D with 6 7^ 0, there exist elements q,r E D such that 

a = bq + r, (36) 

where either r = or v{r) < v{b). This means that in an Euclidean domain division 
with remainder is possible. An example for an Euclidean domain is given by the integer 
numbers Z. 

Finally, a field is a commutative ring in which every nonzero element has a multiplica- 
tive inverse, e.g. R — {0} is an abelian group. Any field is automatically an Euclidean 
domain. Examples for fields are given by the rational numbers Q, the real numbers M, the 
complex numbers C or Zp, the integers modulo p with p a prime number. 
Fig. (5) summarizes the relationships between the various domains. Of particular impor- 
tance are polynomial rings in one or several variables. Fig. (6) summarizes the structure 
of these domains. Note that a multivariate polynomial ring R[ always be 

viewed as an univariate polynomial ring in one variable Xn with coefficients in the ring 

R[xi, ...,Xn-l]. 
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The algorithm for the calculation of the gcd in an Euclidean domain dates back to Euclid 
[30]. It is based on the fact that if a = 6g + r, then 

gcd(a,6) = gcd(6,r). (37) 

This is easily seen as follows: Let c = gcd(a, b) and d = gcd(6, r). Since r — a — bq we see 
that c divides r, therefore it also divides d. On the other hand d divides a — bq + r and 
therefore it also divides c. We now have c\d and d\c and therefore c and d are associates. 
It is clear that for r = 0, e.g. a = bq we have gcd(a, b) = b. Let us denote the remainder 
as r = rem(a, b). We can now define a sequence = a, Vi = b and rj = rem(rj_2, Tj-i) for 
i >2. Then there is a finite index k sucht that r^+i — (since the valuation map applied 
to the remainders is a strictly decreasing function) . We have 

gcd(a, b) = gcd(ro, n) = gcd(ri, r2) = ... = gcd(rfc-i, r^) = r^. (38) 

This is the Euchdean algorithm. We briefly mention that as a side product one can find 
elements s, t such that 

sa + tb — gcd(a, 6). (39) 
This allows the solution of the Diophantine equation 

sa + tb = c, (40) 
for s and t whenever gcd(a, b) divides c. 

We are primarily interested in gcd computations in polynomial rings. However, polyno- 
mial rings are usually only unique factorization domains, but not Euclidean domains, e.g. 
division with remainder is in general not possible. As an example consider the polynomials 

a{x) = + 2x + 3 and b{x) = 5a; + 7 in 7j[x]. It is not possible to write a{x) in the 
form a{x) = b{x)q{x) +r(a;), where the polynomials q{x) and r{x) have integer coefficients. 
However in 'Q[x\ we have 

x^ + 2. + 3 = (5x + 7)(lx+l)+| (41) 

and we see that the obstruction arrises from the leading coefficent of b{x). It is therefore 
appropriate to introduce a pseudo-division with remainder. Let D[x] be a polynomial ring 
over a unique factorization doamin D. For a{x) = a.^x'^ + ... + Oq, b{x) — b^x'^ + ■■■ + &o 
with n>m and b{x) ^ there exists q{x),r{x) e D[x] such that 

6r"^+^a(x) = b(x)q(x)+r(x) (42) 

with deg(r(a;)) < deg{b{x)). This pseudo-divison property is sufficient to extend the Eu- 
clidean algorithm to polynomial rings over unique factorization domains. 
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Unfortunately, the Euclidean algorithm as well as the extended algorithm with pseudo- 
division have a severe drawback: Intermediate expressions can become quite long. This 
can be seen in the following example, where we would like to calculate the gcd of the 
polynomials 

a{x) = x^ + x^ -3x^-3x^ + Sx'^ + 2x-5, 

b{x) = 3a;^ + 5x^-4x^-92; + 21, (43) 

in Z[x]. Calculating the pseudo- remainder sequence ri{x) we obtain 

r^ix) = -15x^ + 3x^-9, 

rsix) = 15795x^ + 30375X - 59535, 

uix) = 1254542875143750X - 1654608338437500, 

r5(x) = 12593338795500743100931141992187500. (44) 

This implies that a{x) and b{x) arc relatively prime, but the numbers which occur in the 
calculation are large. An analysis of the problem shows, that the large numbers can be 
avoided if each polynomial is split into a content part and a primitive part. The content 
of a polynomial is the gcd of all it's coefficients. For example we have 

15795x2 + 30375X - 59535 = 1215 (l3x2 + 25x + 49) (45) 

and 1215 is the content and 13x^ + 25x + 49 the primitive part. Taking out the content of 
a polynomial in each step requires a gcd calculation in the coefficient domain and avoids 
large intermediate expressions in the example above. However the extra cost for the gcd 
calculation in the coefficient domain is prohibitive for multivariate polynomials. The art 
of gcd calculations consists in finding an algorithm which keeps intermediate expressions 
at reasonable size and which at the same time does not involve too much computational 
overhead. An acceptable algorithm is given by the subresultant method [31, 32]: Similar 
to the methods discussed above, one calculates a polynomial remainder sequence ro(x), 
ri(x), ... rfc(x). This sequence is obtained through ro(x) = a(x), ri(x) = b{x) and 



cf+Vj_i(x) = qi{x)ri{x) + diri+i{x) (46) 

where q is the leading coefficient of rj(x), 5i — deg(ri_i(x)) — deg(ri(x)) and di — {—lY^~^^, 
di — — ioT 2 < i < k. The V'i are defined hy ijji — —1 and 

= {-c,.,f-^ (47) 



Then the primitive part of the last non- vanishing remainder equals the primitive part of 
gcd(a(x),6(x)). 

I finally discuss an heuristic algorithm for the calculation of polynomial gcds [33]. In 
general an heuristic algorithm maps a problem to a simpler problem, solves the simpler 
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problem and tries to reconstruct the solution of the original problem from the solution of 
the simpler problem. For the calculation of polynomial gcds one evaluates the polynomials 
at a specific point and one considers the gcd of the results in the coefficient domain. Since 
gcd calculations in the coefficient domain are cheaper, this can lead to a sizeable speed-up, 
if both the evaluation of the polynomial and the reconstruction of the polynomial gcd can 
be done at reasonable cost. Let us consider the polynomials 

a{x) = + 21a;3 + 35a;^ + 27x + 7, 

b{x) = 12x^-3x^-17x^-45x + 21. (48) 

Evaluating these polynomials at the point ^ = 100 yields a(lOO) = 621352707 and 6(100) — 
1196825521. The gcd of theses two numbers is 

c = gcd(621352707, 1196825521) = 30607. (49) 

To reconstruct the polynomial gcd one writes c in the ^-adic representation 

c = co + Cie + ... + c„e", -|<Ci<|- (50) 
Then the candidate for the polynomial gcd is 

g{x) = C0 + C1X + ... + Cnx"'. (51) 

In our example we have 

30607 = 7 + 6- 100 + 3- 100^ (52) 

and the candidate for the polynomial gcd is g{x) = 3x'^ + 6a; + 7. A theorem guarantees 
now if ^ is chosen such that 

e>l + 2min(||a(x)||oo,||K^)||oo), (53) 

then g{x) is the greatest common divisor of a{x) and b{x) if and only if g{x) divides 
a{x) and b{x). This can easily be checked by a trial division. In the example above, 
g{x) = 3x^ + 6x + 7 divides both a{x) and b{x) and is therefore the gcd of the two 
polynomials. 

Note that there is no guarantee that the heuristic algorithm will succeed in finding the gcd. 
But if it does, this algorithm is usually faster than the subresultant algorithm discussed 

previously. Therefore the strategy employed in computer algebra systems like Maple or 
GiNaC is to try first a few times the heuristic algorithm with various evaluation points and 
to fall back onto the subresultant algorithm, if the gcd has not been found by the heuristic 
algorithm. 
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original problem 



solution of original problem 



reduction 



reconstruction 



computation 



simpler problem 



solution of simpler problem 



Figure 7: The modular approach: Starting from the original problem, one first tries to find 
a related simpler problem. The solution of the simpler problem is used to reconstruct a 
solution of the original problem. 

3.4 Remarks 

In this section we discussed issues related to efficiency. It is worth to recall two general 
strategies: 

• The first one is the "divide-and-conquer" approach. If a large task can be partitioned 
recursively into smaller units, this can lead to a considerable speed-up. We have seen 
examples for this approach in Karatsuba's algorithm, in the algorithm for the fast 
Fourier transform and in the search for elements in a binary tree. 

• The second strategy is the modular technqiue. Instead of solving a complicated 
problem right from the start, one tries first to find a related simpler problem. From 
the solution of the simpler problem, one tries then to reconstruct the solution of the 
original problem. Since the reduction to the simpler problem might involve some 
information loss, the reconstruction of the full solution may not be possible or may 
not be unique. However, if there is a simple and efficient way to check, if a guess for 
a solution is the correct solution, this method can be highly competitive. Fig. (7) 
summarizes this approach in a diagram. We discussed as an example for this strategy 
the heuristic gcd algorithm. 



The development of computer algebra systems has triggered research in mathematics on 
constructive algorithms for the solution of some important problems. Examples are the 
factorization of polynomials, symbolic integration and symbolic summation, and simplifica- 
tions with the help of Grobner bases. Most of these topics involve a big deal of mathematics. 
Here I can only sketch an outline of the algorithms. 
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Classical Algorithms 
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4.1 Factorization 



It occurs frequently that we would like to factorize a polynomial. Here I discuss the fac- 
torization of a polynomial u{x) G Z[x] over the integers Z. I outline one specific algorithm 
due to Berlekamp. The algorithm for the factorization is divided into three steps. The 
first (and relative simple) step factors out the greatest common divisior of the coefficients 
and performs a square-free decomposition. The second step uses Berlekamp's algorithm 
and factors the polynomial in the ring Zp, where p is a prime number. If p were sufficiently 
large, such that all products of coefficients would always lie between — p/2 and p/2, the 
factorization over Z could be read off from the factorization over Zp. Unfortunately, the 
required bound on p can become rather large. One uses therefore a small prime number 
p and reconstructs in step 3 from a factorization over Zp a factorization over Zpr . Tising a 
technique called Hensel lifting. Since r can be chosen large, this will yield a factorization 
over Z. 

It should be noted that there are also other (and more efficient) algorithms for the factoriza- 
tion of polynomials. Examples are the probabilistic algorithm by Cantor and Zassenhaus 

[34] or the algorithm by Kaltofcn and Shoup [35]. The computer program NTL by Shoup 
[27] provides a state-of-the-art implementation for the factorization of univariate polyno- 
mials. 

4.1.1 Square- free decomposition 

Suppose a polynomial u{x) contains a factor v{x) to some power m, e.g. 

u{x)^iv{x)rr{x). (54) 

The exponent m can easily be reduced to one as follows: One starts by forming the 
derivative 

u'{x) = m {v{x))"''~^ v'{x)r{x) + {v{x))"^ r'{x) (55) 
and calculates the gcd 

g{x) = gcd {u{x) , u' {x)) . (56) 

Prom eq. (54) and eq. (55) one sees that {v{x))^''^ is a factor of the gcd and u{x) can be 
written as 

Note that g{x) divides u{x) by the definition of the gcd. This process can be repeated, 
until in each term all factors occur with power one. The computational cost is one (or 
several) gcd calculation(s). 
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4.1.2 Berlekamp's algorithm 

After square-free decomposition we can assume that our polynomial 

u{x) = pi{x)p2{x)...pk{x) (58) 

is a product of distinct primes. We assume that deg u{x) = n. Berlekamp's algorithm [36] 
factorizes this polynomial in Zp, where p is a prime number. For < A; < n one obtains 
coefficients Qkj from the mod u{x) representation of x'^^, e.g. 

x''^ = {qk,n-ix"'~^ + ... + qk,ix + qk,o) mod u{x). (59) 
This defines an x n matrix 

(Q0,Q QO,! ... Q0,n-l \ 

(60) 
Qn-1,0 Qn~l,l ... Qn-l,n-l ) 

A non-trivial solution of 

{vQ,vx,...,Vn-\)Q = ('^;o,^;l, ...,^;n-l) (61) 
defines then a polynomial v{x) by 

v{x) = Vn_iX'^~^ + ... + ViX + Vq. (62) 

Note that the solution of eq. (61) is obtained by linear algebra in Zp. The calculation of 

gcd{u{x),v{x) — s) (63) 
for < s < p will then detect the non-trivial factors of u{x) in Zp. 

4.1.3 Hensel lifting 

Assume now that we have a factorization 

u{x) — Vi{x)wi{x) mod p. (64) 

The Hensel lifting promotes this factorization to a factorization 

u{x) — Vr{x)wr{x) mod . (65) 

Since u{x) was assumed to be square-free, we have gcd(f i(a;), ^1(2;)) = 1 mod p. For 
simplicity we also assume that the leading coefficient of u{x) is 1. 

Since gcd(f i(a;), ■u;i(a:)) = 1 mod p we can compute with the Euclidean algorithm polyno- 
mials a{x) and b{x) with deg a{x) < deg wi{x) and deg b{x) < deg vi{x) such that 

a{x)vi{x) + b{x)wi{x) — Imodp. (66) 



31 



Suppose now that we arc given {vr,Wr) and that we wish to compute {vr+i,Wr+i). To this 
aim one computes first a polynomial Cr{x) such that 

p^Cr{x) = Vr{x)wr{x) — u{x) mod p^~^^ . (67) 

By polynomial division of a{x)cr{x) by Wi{x) one obtains a quotient qr{x) and remainder 
ar{x): 

a{x)cr{x) — qr{x)wi{x) + ar{x) mod p. (68) 

One further sets 

br{x) = b{x)cr{x) + qr{x)vi{x) mod p. (69) 
Vr+i and Wr+i are now given by 

Vr+l — Vr{x) —p^hr{x) modp''"^''^, 

Wr+i — Wr{x) — p^ ar{x) mod p^'^^ . (70) 

As an example let us consider u{x) — x"^ -\- Tlx + 176. This polynomial is already square- 
free and using Berlekamp's algorithm we try to find a factorization in Z3. We use the 
asymmetric representation {0,1,2} for Z3 instead of the more symmetric representation 
{— 1, 0, 1}. We find the factorization 

u{x) — x^ -\-2 mod 3 

= (x + l)(x + 2) mod3. (71) 

Hensel lifting yields then 

u{x) = (x + 7)(a; + 2) mod 9 

= (a; + 16)(a; + ll) mod 27. (72) 

The factorization in Z27 agrees already with the factorization in Z. 
4.2 Symbolic integration 

Textbook integration techniques are often not more than a collection of heuristic tricks, 
together with a look-up integration table. A major achievement, which was initiated by 
the development of computer algebra systems, was the invention of a systematic procedure 
to decide whether a given function from a certain class of functions has an integral inside 
this class. If this is the case, the procedure provides also a constructive method to calculate 
this integral. This is the Risch integration algorithm [37]. It is beyond the scope of these 
notes to describe this algorithm in all details, but some important features can already be 
seen, if we restrict ourselves to rational functions. 

A rational function is a quotient of polynomials in the integration variable x: 

_ a{x) _ amx"^ + ... + aix + ap 
■'^^^ ~ h{x) ~ hnx"" + ... + bix + bo ' ^ ^ 
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Suppose that we know the factorization of the denominator over C: 



6„x" + ... + M + &o = bnix - Ci)'^\..{x - Cr^^ (74) 

with mi + ... + mr = n. Then after polynomial division with remainder and partial fraction 
decomposition we can write f{x) in the form 

r mi , 

fix) = pW + ££ (76) 

where p{x) is a polynomial in x and the dij are complex numbers. The integration of p{x) 
is trivial. For the pole terms we have 



dx ( In (x — q) , J = 1, 

{x - CiY \ (lij) {x-ly-^ ' 3 > 1- 



(76) 



The major inconvenience for this approach is the need to factor the denominator completely 
and thereby introducing algebraic extensions (like square roots or complex numbers), even 
if the integrand and the integral can entirely be expressed in a smaller domain (like Q). 
As an example consider the integral 



^2 ^ 



[1 + xY 1 + X 

The method above rewrites the integrand as 

1-x^ 1 



(1+^2)2 2{x + iy 2(x-i)2' 



(78) 



and introduces the complex unit i, which cancels out in the final result. 
As can be seen from eq. (75) and eq. (76) the result for the integration of a rational 
function can always be written as a rational function plus some logarithms. It is desirable 
to compute as much of the integral as possible in the domain of the integrand, and to find 
the minimal algebraic extension necesary to express the integral. This is done in two step: 
The first step uses Hermite's reduction method and reduces the problem to one where the 
denominator is square-free. The second step apphes then the Rothstein-Trager algorithm 
to obtain the logarithmic part. We first consider Hermite's reduction method. Let 

fix) = P(^) + ^, (79) 

b{x) 

with deg a{x) < deg b{x) and gcd(a(a;), 6(a;)) = 1. Assume further that m is the highest 
power to which irreducible factors occur in the factorization of the denominator, e.g. 

b{x) = u{x){v{x)r. (80) 
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With the help of the Euchdean algorithm we may then calculate polynomials r{x) and s{x) 
such that 

r{x)u{x)v'{x) + s{x)v{x) — — a{x). (81) 

Then we obtain for the integral 

f _ rjx) f ^ Jl-m)s{x)-u{x)r'{x) 

J ''\{x){v{x)r ~ {v{x)r-' J u{x){v{x)r-' ' ^ ^ 

and the problem is reduced to one with a smaller power of v{x) in the denominator. 
Therefore wc can assume that wc arc left with an integral of the form f[x) = a{x)/b{x), 
where deg a{x) < deg b{x) and b{x) is square-free. This integral will yield a result of the 
form 

dx^^ = In (x - Ci) (83) 

where the the zeros of b{x) and the rj's are the residues of f{x) at the q's. There 

are efficient algorithms for the determination of the constants q and rj. The first algorithm 
was invented by Rothstein [38] and Trager [39]. The method was later improved by Trager 
and Lazard and Rioboo [40]. 

Before concluding this section, let us comment on the Risch integration algorithm. The 
integration method outlined above for the integration of rational functions was generalized 
by Risch to an algorithm for elementary functions. The class of elementary functions is de- 
fined as follows: An elementary function of a variable x is a function that can be obtained 
from rational functions by repeatedly adding a finite number of logarithms, exponentials 
and algebraic numbers or functions (for example algebraic numbers like the complex unit i 
or square roots). Any (finite) nested combination of those functions is again an elementary 
function. 

For a systematic approach to integration it is useful to reformulate some concepts in an al- 
gebraic way: A differential field is a field F of characteristic with a mapping D : F ^ F, 
which satisfies 

D{f + g) = D{f)+D{g), 
D(f-g) = D(f)-g + f.D(g). (84) 

The mapping D is called a derivation or differential operator. An example for a differen- 
tial field is given by the field of rational functions in one variable with D being the usual 
differentiation. 

We have seen in example eq. (77) that the integral for this example is again a rational 
function. In general, however, the integral of a rational function is expressed by logarithms 
and rational functions. In this case, the logarithms are an extension to the original differ- 
ential field of rational functions. In the case of elementary functions, we have to consider 
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logarithmic, exponential and algebraic extensions. 

Let F be a differential field and let G be a differential extension of it. A function g e G is 
called logarithmic over F, if there exists a / e F such that 

9' = f (85) 
A function gf e G is called exponential over F, if there exists a / e F such that 

^ = /'. (86) 

A function g e G is called algebraic over F, if there exists a polynomial p e F[z] such that 

Pig) = 0. (87) 

Given an elementary function /, the Risch algorithm will decide whether the integral of 

/ can be expressed as an elementary function. If this is the case, the algorithm provides 
also a constructive method to express the integral in terms of elementary functions. If 
this is not the case, we know at least, that the integral cannot be expressed in terms of 
elementary functions. 



4.3 Symbolic summation 

Symbolic summation can be thought of as the discrete analog of symbolic integration. For 
indefinite integration one looks for a function F{x) such that F'{x) = f{x). Then the 
definite integral over f{x) is given by 

b 

J dxf{x) = F(6) -F(a). (88) 

a 

For indefinite summation one looks for a function T{k) defined for integers k, such that 
the difference 



T{k + 1)-T{k) = t{k) (89) 
equals a given function t{k). The definite sum is then given by 

6-1 

= T{b)-T{a) (90) 

k=a 

One now looks for an algorithm, which for a given class of functions f{k) decides, if the 
indefinite sum can be expressed in terms of functions of this class. If the answer is positive. 
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we would like to have a function T{k), which fuUfills cq. (89). This problem has been solved 
for the class of "hypergeometric terms" . A sum of hypcrgeometric terms 

n-l 

J2m (91) 

k=0 

is a sum for which the ratio of two consecutive terms is a rational function r{k) of the 
summation index k, e.g. 

A function t{k), which satisfies eq. (92) is called a hypergeometric term. In the following 
we will outline Gosper's algorithm [41]. Suppose first that we can write r{k) in the form 

where a{k), b{k) and c{k) are polynomials in k and 

gcd{a{k),b{k + j)) = 1, (94) 

for all non- negative integers j. In fact, such a factorization exists for every rational function 
r{k) and there is a systematic way to find it. We then look for a rational function x{k), 
which satisfies the first order difference equation 

a(k)x{k + 1) - b(k - l)x(k) = c(k). (95) 

If such a solution x{k) does not exist, the sum in eq. (91) cannot be done within the class 
of hypergeometric terms. Otherwise, the solution is given by 



n—l 



E,, , b(n — l)x(n) , , 
m = ^ . ^ ^ ^w- (96) 

fc=0 ^ ' 

This completes Gosper's algorithm. Up to now we considered indefinite summation. Al- 
though a solution for an indefinite summation problem may not exist, we might be able 
to find a solution for a specific definite summation problem. Definite summation problems 
have been studied for sums of the type 

oo 

/(n) = ^(^'^)' (9^) 

fe=— oo 

where both 

F(n + 1.A;) , F{n,k^\) 

—— — — ^ and — — ^ (98) 

F{n,k) F{n,k) ^ ' 

arc rational functions of n and k. Algorithms for these type of sums have been given 
by Sister Celine, Wilf and Zeilberger and Pctkovsck. The book by Petkovsek, Wilf and 
Zeilberger [42] gives a good introduction to this subject. 
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4.4 Grobner bases 



In this paragraph we consider simpHfications with the help of Grdbner bases. Assume that 
we have a (possibly rather long) expression /, which is a polynomial in several variables 
Xi, Xfc. In addition we have several siderelations of the form 

Sj{xi,...,Xk) = 0, l<j<r, (99) 

which are also polynomials in xi, Xk. A standard task is now to simplify / with respect 
to the siderelations Sj, e.g. to rewrite / in the form 

/ = GiSi + ... + QrSr + g, (100) 

where g is "simpler" than / The precise meaning of "simpler" requires the introduction 
of an order relation on the multivariate polynomials. As an example let us consider the 
expressions 

h = x + 2y\ f2 = x\ (101) 
which we would like to simplify with respect to the siderelations 

si = x^ -\- 2xy, 

52 = xy + 2y^-l. (102) 

As an order relation we choose lexicographic ordering, e.g. x is "more complicated" as y, 
and x'^ is "more complicated" than x. This definition will be made more precise below. 
A naive approach would now take each siderelation, determine its "most complicated" 
element, and replace each occurence of this element in the expression / by the more simpler 
terms of the siderelation. As an example let us consider for this approach the simplification 
of /2 with respect to the siderelations Si and 

/2 = = si - 2xy = si - 2ys2 + V - 2?/, (103) 

and /2 would simplify to 4|/'^ — 2y. In addition, since /i does not contain x"^ nor xy, the 
naive approach would not simplify /i at all. However, this is not the complete story, since 
if si and S2 are siderelations, any linear combination of those is again a valid siderelation. 
In particular, 

53 = ysi — XS2 — X (104) 

is a siderelation which can be deduced from Si and 82- This implies that /2 simplies to 
with respect to the siderelations Si and S2. Clearly, some systematic approach is needed. 
The appropriate tools are ideals in rings, and Grobner bases for these ideals. 
We consider multivariate polynomials in the ring R[xi, x^]. Each element can be written 
as a sum of monomials of the form 

c<^..a;^^ (105) 
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We define a lexicographic order of these terms by 



cx 



mi 



.X 



> cx 



■x. 



(106) 



if the leftmost nonzero entry in (mi — m'-^, 
can write any element / e Xfe] as 



., ruk — m'lJ is positive. With this ordering we 



(107) 



i=0 



where the hi are monomials and /ij+i > hi with respect to the lexicographic order. The term 
hn is called the leading term and denoted lt(/) = hn- Let B = {bi, ...,br} C R[xi, ...,Xk] 
be a (finite) set of polynomials. The set 



^ttik 



tti e R[xi, ...,Xk] 



(108) 



is called the ideal generated by the set B. The set B is also called a basis for this ideal. 
(In general, given a ring R and a subset 7 C -R, / is called an ideal if a + & e 7 for all 
a,b,^ I and ra E I for all a G / and r E R. Note the condition for the multiplication: The 
multiplication has to be closed with respect to elements from R and not just I.) 
Suppose that we have an ideal / and a finite subset H <Z I. We denote by It (if) the set 
of leading terms of H and, correspondingly by lt(/) the set of leading terms of I. Now 
suppose that the ideal generated by lt{H) is identical with the one generated by lt(/), e.g. 
It (if) is a basis for (lt(/)). Then a mathematical theorem guarantees that H is also a basis 
for /, e.g. 



{lt{H)) = (lt(7)) 



{H)^I 



(109) 



However, the converse is in general not true, e.g. if ii is a basis for / this docs not imply 
that lt{H) is a basis for (lt(/)). A further theorem (due to Hilbert) states however that 
there exists a subset G C I such that 



{G) = I and (It(G')) = (lt(/)). 



(110) 



e.g. G is a basis for / and It(G') is a basis for (lt(/)). Such a set G is called a Grobncr basis 
for /. Buchberger [43] gave an algorithm to compute G. The importance of Grobner bases 
for simphfications stems from the following theorem: Let G be a Grobner basis for an ideal 
I C R[xi, ...,Xk\ and / e R[xi, ...,Xk\- Then there is a unique polynomial g e R[xi, ...,Xk\ 
with 



;iii) 



and no term of g is divisible by any monomial in lt(G). 

In plain text: / is an expression which we would like to simplify according to the siderela- 
tions defined by I. This ideal is originally given by a set of polynomials {si, Sr} and the 



38 



siderelations are supposed to be of the form Si = 0. From this set of siderelations a Grobner 
basis {61, ...,br'} for this ideal is calculated. This is the natural basis for simplifying the 
expression /. The result is the expression g, from which the "most complicated" terms of 
G have been eliminated, e.g. the terms It(G'). The precise meaning of "most complicated" 
terms depends on the definition of the order relation. 

In our example, {si, S2} is not a Grobner basis for (si, S2), since It(si) = x"^ and lt(s2) = xy 
and 

\t(ysi-xs2)^x ^ (lt(si),lt(s2)). (112) 
A Grobner basis for (si, S2) is given by 

{x,2y'-l}. (113) 
With bi — X and 62 = — 1 as a Grobner basis, /i and /2 can be simplified as follows: 

/i = 61 + 62 + 1, 

/2 = a;6i + 0, (114) 
e.g. /i simplifies to 1 and /2 simplifies to 0. 



4.5 Remarks 

The Risch integration algorithm and Gosper's algorithm for summation both operate on a 
class of functions and can decide if the integration- or summation problem can be solved 
inside this class of functions. If this is the case, they also provide the solution in terms of 
functions of this class. The relevant classes of functions are the class of elementary functions 
for integration and the class of hypergeometric terms for summation. Unfortunately, these 
classes do not contain important functions relevant to loop calculations in perturbative 
quantum field theory. For example, Euler's dilogarithm 

U2{x) ^ -J dt^^^^^^ (115) 



is not in the class of elementary functions. For the class of hypergeometric terms, the 
situation is even worse. The harmonic numbers 

i=i 

are not included in the class of hypergeometric terms. Therefore, logarithms like 

00 j 

-ln(l-x) ^ 
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oo 



U2{x) 



■2 ' 



^ln^(l-x) 



oo j i—1 



J 



1 



(117) 



1=1 j=i 



are not in the class of hyper geometric terms. The three functions in eq. (117) are the three 
basic transcendental functions occuring in one-loop calculations. Beyond one loop, there 
are additional functions. 



In this section I discuss algorithms relevant to perturbative calculations in high energy 
physics. Topics include: Generation of all contributing Feynman diagrams, contraction 
of indices and Dirac algebra gymnastic, reduction of tensor loop integrals to scalar loop 
integrals and the evaluation of scalar loop integrals in terms of analytic functions. This 
section can only offer a selection of topics. Not included are for example methods based on 
Mellin transformations [44] or asymptotic expansions. The last method has been reviewed 
in [45, 46]. 

5.1 Graph generation 

The first step in any perturbative calculation is usually the determination of all relevant 
Feynman diagrams. This is trivial for processes involving not more than a handful dia- 
grams, but requires a systematic procedure for processes involving a few hundreds or more 
diagrams. An efficient algorithm has been developed by Nogueira [47] and implemented 
into the program QGRAF [48]. The algorithm avoids recursive generation of Feynman di- 
agrams and comparisons of stored diagrams. The algorithm for the generation of Feynman 
graphs is divided into two steps: In the first step all relevant topologies are generated. A 
topology is just a collection of nodes and edges, which connect nodes. A topology with 
n nodes is represented by a n x n adjacency matrix A. The entry Aij of the adjacency 
matrix denotes the number of edges connecting the nodes i and j. The algorithm generates 
all relevant adjacency matrices. Using an order relation for adjacency matrices, one can 
avoid to generate similar adjacency matrices, e.g. matrices which may be related to each 
other by a permutation of the nodes. In the second step the external nodes are labelled 
with the external particles and the edges and internal nodes are substituted by propagators 
and interaction vertices in all possible ways compatible with the Feynman rules. To avoid 
to generate equivalent Feynman graphs more than once during this stage, the symmetry 
group of each Feynman graph is determined and a further order relation is used to return 
only the graph which is the "smallest" within its symmetry class with respect to this order 
relation. In addition the symmetry factor is returned. 

Other algorithms for the generation of Feynman diagrams have been considered for exam- 
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pie by Kiiblbeck, Bohm and Denner [49] and implemented into FeynArts [50] as well as by 
Kajantie, Laine and Schroder [51]. 



5.2 Contraction of indices and Dirac algebra 

After all diagrams have been generated, they are translated to mathematical expressions 
with the help of Feynman rules. Edges correspond to propagators and vertices to interac- 
tion vertices. As an example for Feynman rules we give the rule for gluon propagator in a 
covariant gauge and the rule for the quark-gluon vertex: 

^^^^ = 

N (118) 

The resulting expressions involve summations over repeated indices. The contraction of 
these indices can be done by applying succesivly a few rules: 

Pt^q^ = pq, P^l^ = A 7m7^ = D. (119) 
If a contraction over a string over Dirac matrices occurs, like in 

7m7.7^ (120) 
one uses first the anti-commuation relations of the Dirac matrices 



{7^,7,} = 2g^, (121) 

to bring the two matrices with the same index next to each other and uses then eq. (119). 
For the calculation of squared amplitudes there will be always a trace over the strings of 
Dirac matrices, which can be evaluated according to 

Tr 1 = 4, (122) 

Tr 7/Jl7/J2---7At2n ~ 9/J.lfl2'~^^ 7M3"-7At2n ~ S'miMs'^^ 7/J2 7/14---7At2n "I" ••• "I" 9/11/12,1^^ 7/J2 • • •7At2n-l • 

Traces over an odd number of Dirac matrices vansih. The recursive procedure for the 
evaluation of a trace over 2n Dirac matrices will generate (2n — l)(2n — 3)...3-l = (2?7, — 1)!! 
terms. This number grows exponentially with n. If there are no further relations between 
the indices /ii, /i2n this is indeed the number of terms in the final result. However, in 
almost all practical applications, the free indices /zi, /i2n get contracted with a tensor, 
which is at least symmetric in some of its indices. In that case the recursive procedure is 
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inefficient. A better way consists in splitting the string of Dirac matrices into two parts 
[52]: 

) (7^,+i-7^2j (123) 

Tfiis process is repeated until each factor contains only a few Dirac matrices (say not more 
than two) . One then rewrites each factor in terms of the basis 

r(o) = 1, 

rll = ^[7.1,7.2], etc. (124) 
The general term of this basis has the form 

^uL-i^, = lluilu2--lu,], (125) 

where [...] denotes anti-symmetrization. Obviously, 

T^^l'L.u, = forj>l. (126) 
Products of these basis elements can be combined with a Clebsch-Gordan type formula: 

pW _ npi-Pk pW ('1271 

k=\i-j\ 

The Clebsch-Gordan coefficients C^j"'^^j,^ ^. are anti-symmetric in the sets {fii...fii}, 
and {pi-.-Pk}- Consider now the situation where the tensor with which the free indices /xi, 
pL2n s-re contracted, is symmetric in /ij and pj. In this situation, any term involving 
a Clebsch-Gordan coefficient with fii and pj in the same index field can immediately be 
discarded. This can lead to a considerable speed-up. 



Minor complications occur if 75 appears in the calculation. Within dimensional regu- 
larisation this requires a consistent definition of 75. One possible definition is the 't Hooft - 
Veltman scheme, which takes 75 as a generic four-dimensional object. As a consequence one 
has to distinguish between four- dimensional and D = A — 2e dimensional quantities. In this 
scheme D is assumed to be greater than 4. It is further assumed that the four-dimensional 
subspace and the (— 2^) dimensional subspace are orthogonal. For the calculation one splits 
all quantities into a four-dimensional part (denoted with a tilde) and a (—2s) dimensional 
part (denoted with a hat) . 

9,^u = 9^,u + g„i., 7m = 7m + 7m- (128) 
These quantities satisfy relations like 

^m''^^, ^/ = -2£, ~g,.r'-0. (129) 
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75 is then defined as a generic four- dimensional object: 

75 = jfa^^sTfl'l'- (130) 

As a consequence, 75 anti-commutes with the four-dimensional Dirac matrices, but com- 
mutes with the remaining ones: 

{75,%} = 0, [75,7/.] =0. (131) 

The program FORM [13] is one of the most efficient tools for manipulations involving 
contraction of indices and traces over Dirac matrices. The program TRACER [53] has 
been written for the manipulations of the Dirac algebra in the 't Hooft - Veltman scheme. 



5.3 One-loop integrals and Passarino- Veltman reduction 

We now consider the reduction of tensor loop integrals (e.g. integrals, where the loop 
momentum appears in the numerator) to a set of scalar loop integrals (e.g. integrals, 

where the numerator is independent of the loop momentum). For one- loop integrals a 
systematic algorithm has been first worked out by Passarino and Veltman [54]. Consider 
the following three-point integral 

where pi and p2 denote the external momenta. The reduction technique according to 
Passarino and Veltman consists in writing I^^ in the most general form in terms of form 
factors times external momenta and/or the metric tensor. In our example above we would 
write 

= PKC21+P^P^C22 + K,P^}C23 + ^'^''C24, (133) 

where {^1,^2} = P1P2 +P2Pi- One then solves for the form factors C2i,C22,C23 and C24 
by first contracting both sides with the external momenta PiPi, ^2^2, {PiiVVi 
metric tensor gi^^ . On the left-hand side the resulting scalar products between the loop 
momentum and the external momenta are rewritten in terms of the propagators, as for 
example 

2pi • A; = k'^-{k- pif + pi- (134) 

The first two terms of the right-hand side above cancel propagators, whereas the last 
term does not involve the loop momentum anymore. The remaining step is to solve for 
the formfactors €2% by inverting the matrix which one obtains on the right-hand side of 
equation (133). Due to this step Gram determinants usually appear in the denominator of 
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Pi. 



\ 



Figure 8: An example for irreducible scalar products in the numerator: The scalar product 
2pik2 cannot be expressed in terms of inverse propagators. 



the final expression, 
the triangle 



In the example above we would encounter the Gram determinant of 



A, 



pI Pi ■ P2 
Pi ■ P2 pI 



(135) 



One drawback of this algorithm is closely related to these determinants : In a phase space 
region where pi becomes coUinear to p2, the Gram determinant will tend to zero, and the 
form factors will take large values, with possible large cancellations among them. This 
makes it difficult to set up a stable numerical program for automated evaluation of tensor 
loop integrals. There are modifications of the Passarino-Veltman algorithm, which avoid 
to a large extent the appearance of Gram determinants. These improved algorithms are 
based on spinor methods [55]. 



5.4 Beyond one- loop 

The Passarino-Veltman algorithm is based on the observation, that for one-loop integrals 
a scalar product of the loop momentum with an external momentum can be expressed as 
a combination of inverse propagators. This property does no longer hold if one goes to 
two or more loops. Fig. (8) shows a two-loop diagram, for which the scalar product of 
a loop momentum with an external momentum cannot be expressed in terms of inverse 
propagators. 

A different and more general method for the reduction of tensor integrals is needed. In 
the following I will outline algorithms, which work beyond one loop. The first step is 
to relate tensor integrals to scalar integrals with raised powers of the propagators and 
different values of the space-time dimension D [56]. The starting point is the Schwinger 
parameterization for the propagators 

oo 

(Z^ = ^/rfxx^-^exp(xP). (136) 



Let us first consider a scalar two-loop integral. Combining the exponentials arising from 
different propagators one obtains a quadratic form in the loop momenta. For instance, for 
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a given two- loop integral with loop momenta ki and /c2, one has then 




r d^ki r d^k2 1 



(137) 



and 



n 



a kl + b kl + 2 cki ■ k2 + 2 d ■ ki + 2 e ■ k2 + f . 



(138) 



1=1 



The momenta k^, . . . ,kk are linear combinations of the loop momenta ki, /c2 and the external 
momenta. The coefficients a, b, c, d'^, and / are directly readable from the actual graph: 
a{b) = where the sum runs over the legs in the ki {k2) loop, and c = with 

the sum running over the legs common to both loops. With a suitable change of variables 
for the loop momenta ki,k2, one can diagonalize the quadratic form and the momentum 
integration can be performed as Gaussian integrals over the shifted loop momenta according 
to 



Let us now consider a tensor integral. After the change of variables for the diagonalization 
of the quadratic form, we have a polynomial in the Schwinger parameters and the loop 
momenta ki and k2 in the numerator. Integrals with an odd power of a loop momentum 
in the numerator vanish by symmetry, while integrals with an even power of the loop 
momentum can be related by Lorentz invariance to scalar integrals: 



The generalization to arbitrary higher tensor structures is obvious. In the remaining 
Schwinger parameter integrals, the tensor integrals introduce additional factors of the pa- 
rameters Xi and oil/V. These additional factors can be absorbed into scalar integrals with 
higher powers of propagators and shifted dimensions, by introducing operators i"*", which 
raise the power of propagator i by one, or an operator d"*" that increases the dimension by 
two. 




(139) 




(140) 



oo 




(141) 
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All Schwinger integrals are rewritten in terms of these scalar integrals. Therefore, using 
an intermediate Schwinger parametrization, we have expressed all tensor integrals in terms 
of scalar integrals. The price we paid is that these scalar integrals involve higher powers 
of the propagators and/or have shifted dimensions. Each integral can be specified by its 
topology, its value for the dimension D and a set of indices, denoting the powers of the 
propagators. In general the number of different integrals is quite large. There are now 
two different methods to proceed: The first method (integration by parts and solutions 
of differential equations) observes that there are non-trivial relations among the obtained 
set of integrals and employs an elimination procedure to reduce this set to a small set of 
so-called master integrals. These master integrals are then evaluated explicitly. 
Since some integrals have to be evaluated explicitly anyway, one may ask if there is an 
efficient algorithm which evaluates directly these integrals for arbitrary D and arbitrary 
set of indices. This is the philosophy of the second method (expansion of transcendental 
functions) . 

5.4.1 Integration by parts and solutions of differential equations 

Let us start with the method, which reduces the (large) number of scalar integrals to 
a (small) set of master integrals. There exist non-trivial relations among various scalar 
integrals. Some of these relations can be obtained from integration-by-part identities [57] 
or the invariance of scalar integrals under Lorentz transformations [58] . Integration-by-part 
identities are based on the fact, that the integral of a total derivative is zero: 

Here, k is the loop momentum, the pj's are the external momenta and v can either be a loop 
momentum or an external momentum. Working out the derivative yields a relation among 
several scalar integrals. In a similar way one obtains relations based on Lorentz-invariance. 
A scalar integral is evidently invariant under an infinitessimal Lorentz transformation, 
parametrized as 

p^'^p^ + 5jf = + Se'^p'' with 5e(; = -Se^ . (143) 
This implies that 

/ d d d d \ 

+ • • • + Pn^ P^:^ Hpi, • • • ,Pn) = , (144) 

where I{pi, ■■■,Pn) is a scalar integral with the external momenta pi. Again, working out 
the derivatives yields a relation among several scalar integrals. 

Each relation is linear in the scalar integrals and in principle one could use Gauss elimina- 
tion to reduce the set of scalar integrals to a small set of master integrals. In practice this 
approach is inefficient. An efficient algorithm has been given by Laporta [59]. The starting 
point is to introduce an order relation for scalar integrals. This can be done in several 
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ways, a possible choice is to order the topologies first: A scalar integral corresponding to 
a topology Ti is considered to be "smaller" than an integral with topology T2, if Ti can 
be obtained from T2 by pinching of some propagators. Within each topology, the scalar 
integrals can be ordered according to the powers of the propagators and the dimension of 
space-time. Laporta's algorithm is based on the fact that starting from a specific topology, 
integration-by-part and Lorentz-invariance relations only generate relations involving this 
topology and "smaller" ones. To avoid to substitute a specific identity into a large number 
of other identities, one starts from the "smallest" topology and generates all relevant re- 
lations for this topology. Inside this class, integrals with higher powers of the propagators 
or higher dimensions are then expressed in terms of a few master integrals. These manip- 
ulations involve only a small subset of the complete system of integrals and relations and 
can therefore be done efficiently. Once this topology is completed, one moves on to the 
next topology, until all topologies have been considered. All integrals, which cannot be 
eliminated by this procedure are called master integrals. These master integrals have to 
be evaluated explicitly. 

To evaluate these integrals the procedure used in [58, 60] consists in finding first for each 
master integral a differential equation, which this master integral has to satisfy. The deriva- 
tive is taken with respect to an external scale, or a ratio of two scales. It turns out that 
the resulting differential equations are linear, inhomogeneous first order equations of the 
form 

^T{y) + f{y)T{y) = g{y), (145) 

where y is usually a ratio of two kinematical invariants and T{y) is the master integral 

under consideration. The inhomogeneous term g{y) is usually a combination of simpler 
master integrals. In general, a first order linear inhomogeneous differential equation is 
solved by first considering the corresponding homogeneous equation. One introduces an 
integrating factor 

M{y) ^e^ f^y^'^y, (146) 

such that T'o(y) = 1/M{y) solves the homogenous differential equation {g{y) — 0). This 
yields the general solution of the inhomogenous equation as 

^^y^ ^ (/ ^^y^^^y^'^y + ^) ' ^^^^^ 

where the integration constant C can be adjusted to match the boundary conditions. This 
method yields a master integral in the form of transcendental functions (e.g. for example 
hypergeometric functions), which still have to be expanded in the small parameter e of 
dimensional regularization. Although this can be done systematically and will be discussed 
in the next paragraph, the authors of [58, 60] followed a different road. As an example 
let us consider the calculation for e^e^ — > 3 jets, which involves three scales. This scales 
can be taken to be su, S23 and S123. It is convenient to have only one dimension- full 
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quantity S123 and to introduce two dimensionless quantities xi = S12/S123 and X2 = 523/5123- 
Factoring out a trivial dimension-full normalization factor, one writes down an ansatz for 
the solution of the differential equation as a Laurent expression in e. Each term in this 
Laurent scries is a sum of terms, consisting of basis functions times some unknown (and to 
be determined) coefficients. This ansatz is inserted into the differential equation and the 
unknown coefficients are determined order by order from the differential equation. This 
apporach will succeed if we know in advance the right set of basis functions. The basis 
functions are a subset of multiple polylogarithms. In sec. (3.1) we already defined in eq. 
(12) the functions G{zi, z^, y). We now sligthly enlarge this set and define G(0, 0; y) 
with k zeros for Zi to z^ to be 

G'(0,...,0;|/) = ^(Ini/)'. (148) 

This permits us to allow trailing zeros in the sequence (^i, z^). We can now define 
two subsets of these functions. The first subset are harmonic polylogarithms [61] for 
which y — xi and all Zi are from the set {0, 1}. The second subset are two-dimensional 
harmonic polylogarithms [60] for which y — X2 and all Zi are from the set {0,xi,l — 
Xi,l}. The ansatz consists in taking the harmonic polylogarithms and the two-dimensional 
harmonic polylogarithms as a set of basis functions for the three-scale problem e+e" 
3 jets. The coefficients for harmonic polylogarithms are rational functions in xi and X2, the 
coefficients of the two-dimensional harmomic polylogarithms may in addition also contain 
(one-dimensional) harmonic polylogarithms. 



5.4.2 Expansion of transcendental functions 

In this paragraph we now discuss the second approach. Instead of reducing first a large set 
of scalar integrals to a small set of master integrals, which have to be worked out explicitly, 
we may ask if there is an efficient way to evaluate the scalar integrals directly [62, 63, 9]. 
The basic observation is that many relevant integrals can be expressed as (possibly nested) 
sums involving Gamma-functions. As an example let us discuss the one-loop triangle with 
two external masses defined by 

T:i2(m,i/i,i/2,i^3;2;i) = (-5123)"'"^^^'''"" / mo / Ax^^i / Ax^^a / h^u3 : (149) 



^7^ 



where /c2 = ki — pi — p2, = k,2 — ps and xi = S12/S123. We inserted a prefactor 
(—•§123)""^^^^'^^^^ to make the integral dimensionless and used the short-hand notation 
^ij — + ^3 for sums of indices. Within this approach the integral is needed for arbitrary 
(integer) powers of the propagators and possibly also in shifted dimensions 6 — 2£, 8 — 2s, 
etc.. The space-time dimension is written as D = 2m — 2£, where m is an integer. It is 
not too hard to derive a series representation for this integral, consisting of a combination 
of hyper geometric functions 2-^1- 

r^. , . r(£ - m ;/23)r(l - £ + m - ;/23)r(m - £ - t/13) 

Tri2 m,i/i,i/2,i^3;a;i) = ^, ^-p. ^-p. 7, ^ (150 

r(i/i)r(z/2)r(i/3)r(2m - 2£ - 1/123) 
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X 

n=0 



X 



m—e—V2z 



T{n + i^i)r(n - £ + m - z/2) T{n + i'^)T{n - m + e + 2/123) 



r(n + 1 + m - e - 2/23) r(n + 1 - m + £ + 1/23) 



One observes that the small paramter e of dimensional regularization appears in the 
Gamma-functions. This expression has now to be expanded into a Laurent series in e. 
For this particular example it is possible to convert the hypergeometric functions into an 
integral representation, to expand the integrand and to perform the resulting integrals. 
However, more complicated topologies have a representation as nested sums, for which a 
useful integral representation is not known. Fortunately, the expansion in £ can be done 
systematically at the level of nested sums. Using the formula 

r(x+l) = xV{x), (151) 

all Gamma functions can be synchronized to the form V{n + as). They are then expanded 
using the formula 

T{n + e) = T{l + e)T{n) 

X (1 + eZi{n - 1) + e^Zii{n - 1) + £^Zin(n - 1) + ... + £"-^Zn...i(n - 1)) . 

(152) 

Here Euler-Zagier sums, introduced in sec. (3.1) in cq. (10) appear. Collecting terms for 
each order in we obtain expressions involving products of Euler-Zagier sums with the 
same upper summation limit. Since Euler-Zagier sums form an algebra, the multiplication 
can be done. After some additional simple manipulations (partial fraction decomposition 
and adjusting of summation limits) we obtain terms of the form 



00 



E^^-.-.-.(^-l)- (153) 



n=l 

Since these are exactly the harmonic polylogarithms [61] 



00 



n=l 



we arc already done. For topologies involving more scales, we have to generalize this 
approach. The basic quantities are no longer Euler-Zagier sums, but Z-sums defined by 



Z{n;mi,...,mk;xi,...,Xk) = V • • • r^- (155) 



n>ii>i2>...>ifc>0 



Glearly, for xi — ...Xk — 1 we recover the Euler-Zagier sums. All algorithms work also for 
this generalization, in particular, there is a multiplication formula. Of particular impor- 
tance are the sums to infinity: 

^imk,-,mi{xk,--,xi) = Z{oo;mi,...,mk;xi,...,Xk) (156) 
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These are called multiple polylogarithms. There is a close relation between the functions 
G{zi, ...,Zk;y) and Li mk,...,mi{xk, Let US denote by Gmi,...,mk{^i, Zk',y) the G- 

function, where mi — 1 zeros preceed zi, m2 — 1 zeros are inserted between zi and Z2 etc.. 
For example G2,3,i{zi, Z2, Z3; y) = ^(0, zi, 0, 0, ^2, ^3; y)- Then 

^mi,...,mfc(-^l) ^k'l y) 
Lirnfc,.--i"^i (-^fe) •• •) -^l) 

showing that the iterated integral representation eq. (12) and the representation as nested 
sums eq. (156) describe the same class of functions. 

5.4.3 The cut technique 

We conclude the discussion of techniques for loop calculations with a method based on 
unitarity [64] . I will discuss this method for one-loop amplitudes. It has been applied to a 
two- loop calculation by Bern et al. [4]. 

Within the unitarity based approach one chooses first a basis of integral functions Jj G J-. 
For one-loop calculations in massless QCD a possible set consists of scalar boxes, triangles 
and two-point functions. The loop amplitude A^°°p is written as a linear combination of 
these functions 

j^ioop ^ J^cJi + r. (158) 

i 

The unknown coefficients q are to be determined, r is a rational function in the invariants, 
not proportional to any element of the basis of integral functions. The integral functions 
themselves are combinations of rational functions, logarithms, logarithms squared and 
dilogarithms. The latter three can develop imaginary parts in certain regions of phase 
space, for example 

Knowing the imaginary parts, one can reconstruct uniquely the corresponding integral 
functions. In general there will be imaginary parts corresponding to different channels (e.g. 
to the different possibilities to cut a one- loop diagram into two parts). The imaginary part 
in one channel of a one-loop amplitude can be obtained via unitarity from a phase space 
integral over two tree-level amplitudes. With the help of the Cutkosky rules we have 

Im^'^^P = Im /" .'^''^^ , , , / — ATA%'^. (160) 



= {-l)'Gm„...,m, (-, l) , (157) 

\Xi X1X2 Xi...Xk J 



Im In 



iO 



Im Li2 ( 1 



-iO 
iO) 



-t - iO) 
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one-loop amplitude under consideration, yl^^^ and A^£'^'^ are tree-level ampli- 
tudes appearing on the left and right side of the cut in a given channel. Lifting eq. (160) 
one obtains 



where "cut free pieces" denote contributions which do not develop an imaginary part in 
this particular channel. By evaluating the cut, one determines the coefficients q of the 
integral functions, which have an imaginary part in this channel. Iterating over all possible 
cuts, one finds all coefficients q. 

One advantage of a cut-based calculation is that one starts with tree amplitudes on both 

sides of the cut, which are already sums of Feynman diagrams. Therefore cancellations and 
simplifications, which usually occur between various diagrams, can already be performed 
before we start the calculation of the loop amplitude. 

In general, a cut-based calculation leaves as ambiguity the ration piece r in eq. (158), 
which can not be obtained with this technique. One example for such an ambiguity would 



This term does not have a cut and will therefore not be detected in a cut-based calcula- 
tion. However, Bern, Dixon, Dunbar and Kosower [64] have proven the following power 
counting criterion: If a loop-amplitude has in some gauge a representation, in which all 
n-point loop integrals have at most n — 2 powers of the loop momentum in the numerator 
(with the exception of two-point integrals, which are allowed to have one power of the loop 
momentum in the numerator), then the loop amplitude is uniquely determined by its cuts. 
This does not mean that the amplitude has no cut-free pieces, but rather that all cut-free 
pieces are associated with some integral functions. 

In particular = 4 super symmetric amplitudes satisfy the power-counting criterion above 
and are therefore cut-constructible. 

QCD does in general not satisfy the power-counting theorem and leaves as an ambigu- 
ity the rational function r. In principle one can obtain the rational piece r by calculating 
higher order terms in e within the cut-based method. At one-loop order an arbitrary scale 
fi'^'^ is introduced in order to keep the coupling dimensionless. In a massless theory the fac- 
tor yu^^ is always accompanied by some kinematical invariant s"*^ for dimensional reasons. 
If we write symbolically 




(161) 



be 




(162) 




(163) 



the cut-free pieces co(so//i^) ^ can be detected at order e: 




(164) 
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In practice, it is more efficient to take into account additional constraints, like the factor- 
ization in coUinear limits [65], to determine the rational piece r. 

5.5 Remarks 

Let us summarize the state-of-the-art of systematic algorithms for the computation of loop 
integrals occuring in perturbative quantum field theory. In sect. 4 we discussed Risch's 
algorithm for symbohc integration and Gosper's algorithm for symbohc summation. These 
algorithms operate on a certain class of functions and provide the solution if it exists within 
this class of functions. The two classes of functions arc the elementary functions for Risch's 
algorithm and the hypergeometric terms for Gosper's algorithm. None of these two classes 
is sufficiently large to contain all functions, which occur in the calculation of loop integrals. 
A counter-example has been given in sect. 4.5 with the three transcendental functions 
(dilogarithm, logarithm and logarithm squared) which occur in one-loop calculations. 
Beyond one-loop, additional transcendental functions occur in loop calculations. From 
the explicit two-loop calculations we now know that the class of functions occuring in 
loop calculations, should contain the multiple polylogarithms. There are now systematic 
algorithms, which allow the calculation of specific loop integrals. As examples we discussed 
an algorithm based on differential equations and an algorithm based on the expansion of 
transcendental functions. The former is based on the manipulation of iterated integrals, 
whereas the latter is based on the manipulation of nested sums. Although the algorithms 
are quite different, their output is within the same class of functions. 
At the end of these lectures, let us look towards the future: Consider the class of rational 
functions extended by the multiple polylogarithms. Is there an algorithm, similar to Risch's 
or Gosper's algorithm, which decides, if a given integral or sum can be expressed inside 
this class ? If so, is there a constructive way to compute the answer ? These are open 
questions. Calculations in perturbative quantum field theory would certainly profit from 
a positive answer to these questions. 
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